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Abstract 

One method to compute the price of an arithmetic Asian option in a Levy driven model is based 
^ on the exponential functional of the underlying Levy process: If we know the distribution of the 

exponential functional, we can calculate the price of the Asian option via the inverse Laplace trans- 
form. In this paper we consider pricing Asian options in a model driven by a general meromorphic 
Levy process. We prove that the exponential functional is equal in distribution to an infinite prod- 
uct of indepedent beta random variables, and its Mellin transform can be expressed as an infinite 
^SJ product of gamma functions. We show that these results lead to an efficient algorithm for computing 

the price of the Asian option via the inverse Mellin-Laplace transform, and we compare this method 
. with some other techniques. 
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1 Introduction 



One of the popular models for stock price dynamics in mathematical finance is based on a geometric 
Levy process St = SqC^", where X is a Levy process started from zero and 5*0 is the initial stock price. In 
this model the interest rate r is typically constant and it is assumed that the measure P is risk-neutral, 
which means that the process Ste'"^^ is a martingale under measure P. We are interested in calculating 
the price of a fixed strike Asian option, which in this model would be given by 



C{So,K,T) := e-'^E 



So e^-du-K^ 



(1) 



In order to calculate C{So, K, T) one could use either a Monte-Carlo simulation approach, or a semi- 
analytical technique. We would like to emphasize two of these techniques, both of which have important 
connections with other branches of probability. The first technique is based on the observation that, 
while the process Zt := exp(Xu)dM is not Markovian, the process 



Zt := {x + Zt)e 



-Xt 



does satisfy Markov property. The process Zt is known as a generalized Ornstein-Uhlenbeck process, 
and it is an important and well-studied object, see [22] and the references therein. Since Z is a Markov 
process, after a change of measure we can rewrite the expectation in (1) so that it involves only Zt, 
and this can be computed by solving the backward Kolmogorov equation, which in this case takes the 
form of a partial integro-differential equation (one dimension for time and one for space variable). This 
approach was developed in a more general setting by Vecer and Xu [30] , and it was skillfully implemented 
by Bayraktar and Xing [2] for jump-diffusion processes. 

The second semi-analytical technique is based on the exponential functional of the Levy process X, 
which is defined for g > as 

e(9) 

I^-= f e^-du. (2) 







Here e{q) denotes an exponential random variable (independent of X) with mean 1/g if g > 0, and 
e(g) = -|-oo if g = 0. In the latter case we need to impose an additional condition E[Xi] < in order to 
ensure that the exponential functional Ig is well-defined (see [3, Theorem 1]). 

Exponential functionals are very important and useful objects, not only in mathematical finance but 
also in many other areas of probability theory. They play a role in such fields as self-similar Markov 
processes, random processes in random environment, fragmentation processes and branching processes. 
They are also connected with the generalized Ornstein-Uhlenbeck processes (their distribution appears 
as a stationary measure of the process). See [3] for an overview of this topic. The connection with Asian 
options is the following: The distribution of the exponential functional can be viewed as the Laplace 
transform in t- variable of the distribution of Zt. By inverting this Laplace transform one can obtain 
the distribution of Zt and, therefore, the price of the Asian option. This approach was use in [6] for 
Levy processes with hyper-exponential jumps (see also [23] for general results on Levy processes with 
one-sided jumps), and this will be the approach that we will follow in the present paper. Our goal is 
to identify the distribution of the exponential functional for a large family of underlying Levy processes 
(having such desirable properties, as jumps of infinite activity or infinite variation) and to apply these 
results to pricing Asian options. 
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Let us say a few words on the current state of knowledge of the distribution of the exponential 
functional. This distribution is known explicitly for processes with hyper-exponential jumps [6], and for 
a slightly more general class of processes with jumps of rational transform [16] (the processes in both 
of these families can have jumps of finite activity only). The distribution of the exponential functional 
is also known for a family of hypergeometric processes, which was introduced recently in [19]. This 
family includes processes with jumps of infinite activity or infinite variation, however in this case the 
distribution of Iq is only known for a single value of g > 0, which depends on the parameters of the 
underlying process X. To the best of our knowledge, there are no known examples of Levy processes 
having double-sided jumps of infinite activity or infinite variation, for which the distribution of Iq can 
be identified explicitly for all g > 0. 

In this paper we fill this gap and we find the distribution of the exponential functional Iq for a general 
meromorphic process. Meromorphic processes (which were introduced in [17]) are defined via the density 
of the Levy measure 



All the coefficients a„, a„, p„ and p„ are strictly positive and the sequences {pn}n>i and {p„}n>i must be 
strictly increasing and unbounded. There are several families of meromorphic processes, for which the 
Laplace exponent iIj{z) := lnE[exp(zXi)] can be computed explicitly, such as beta, theta and hypergeo- 
metric processes, see [14, 15, 18, 19]. The family of meromorphic processes is quite large, in particular 
it is dense in the class of processes with completely-monotone Levy measure (see [17]). At the same 
time, this family is a natural generalization of Levy processes with hyper-exponential jumps (we will call 
them hyper- exponential processes from now on), for which the density of the Levy measure is also given 
by (3), with the two infinite series replaced by finite sums. While hyper-exponential processes can have 
only finite activity compound Poisson jumps, meromorphic processes can exhibit a wide range of path 
behavior, including jumps of infinite activity or infinite variation, which makes them good candidates 
for modeling purposes in mathematical finance. Meromorphic processes can be seen as analogues of 
the well-known models (VG, CGMY, KoBoL, NIG), but with richer analytical structure allowing for 
development of efficient semi-analytical numerical algorithms. See [7, 9] for some recent applications of 
these processes in mathematical finance. 

The paper is organized as follows. In section 2 we study infinite products of independent beta random 
variables. We use these results in section 3 to identify the distribution of the exponential functional of a 
general meromorphic process. In section 4 we discuss numerical issues, such as approximating the Mellin 
transform of the exponential functional and computing the price of an Asian option.. 

2 An infinite product of beta random variables 

In this section we will study infinite products of independent beta random variables. These results will 
be used in section 3, in order to describe the distribution of the exponential functional Iq. The results 
could also be of independent interest, as similar finite products of beta random variables appear in other 
areas of probability, for example, in the definition of the Poisson-Dirichlet distribution, see [27]. 
For a, 6 > 0, let B(a,b) denote the beta random variable, having distribution 




(3) 



n>l 



n>l 




< X < 1. 
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With any two unbounded sequences a = {a„}„>i and /3 = {Pn}n>i which satisfy the interlacing property 

< ai < f3i < a2 < (32 < as < ^3 . . . (4) 
we associate an infinite product of independent beta random variables, defined as 

n>l 



X{a,(3) :=n^K 



(5) 



The random variable X{a, (5) is the main object of interest in this section. Our first task is to establish 
that this random variable is well-defined. 

Proposition 1. The infinite product in (5) converges a.s. 

Before proving Proposition 1, let us establish the following simple (but useful) result. 

Lemma 1. Assume that a and fi are two unbounded sequences satisfying (4), and f : is a 

monotone function such that lima;_j.4.oo f{x) = 0. Then 



J2(fm - /K)) 



n>l 



(6) 



Proof. Assume that / is increasing. Then the condition /(+oo) = implies that f{x) < for all x, thus 
for any m G N we have 



< := X](/(/3n) - /K)) < - /K)) = /Km) - /(«i) < -/(«i) = |/( 



«1 . 



n=l 



n=l 



The above inequality and the fact that the sequence is increasing show that the series in (6) converges 
and its sum is bounded by |/(ai)|. The case when / is decreasing can be considered in exactly the same 
way. □ 



Proof of Proposition 1: Considering the logarithm of both sides of (5), we see that we need to establish 
the a.s. convergence of the infinite series 



log(X(a,/3)) = 5^1n (b^^^ 

n>l ^ 



,/3„-a„) 



/3. 



The Mellin transform of a beta random variable is given by 

T{a + h)T{a + s - 1) 



E 



s-1 



, Re(s) > 1 — a. 



(7) 



(8) 



T{a)T{a + h + s-l] 

By differentiating the above identity twice and setting s = 1, we find 

E [ln(5(,,,))] = V^(a) - ^{a + h), Var [ln(5(,,,))] = ^\a) - + 6), 

where ip{z) := r'{z)/r{z) is the digamma function. It is known that f{z) := ln(2;) —iIj{z) is a completely 
monotone function which decreases to zero (see [1, Theorem 1] or formula 8.361.8 in [12]). This implies 
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that the function —f'{z) + 1/z = ip'iz), has the same property. Applying Lemma 1 we conclude that 
both series 



n>l 

^Var 



In Lb, 



'{a„,l3n-a„) 



n>l 



In LB, 



'(a„,l3n-an) 



(3n 

an 



n>l 

^(V^'K)-v^'(/3„)) 



n>l 



converge, therefore Khint chine- Kolmogorov Convergence Theorem implies a.s. convergence of the infinite 
series (7). □ 



Theorem 1. The Mellin transform Ai{s) := E[(X(a,/3))* exists for all s > 1 — ai and it can be 
analytically continued to a meromorphic function 



M{s) = n 



r(/3„)rK + s - 1) //3„ 



n> 



)r(/3„ + s-i)\ 



(9) 



The above infinite product converges uniformly on compact subsets of the complex plane, not containing 
the poles of Ai{s). 

Proof. Proposition 1 combined with (8) and Levy's Continuity Theorem imply that (9) is true for all 
s G C on the vertical line Re(s) = 1, and that the infinite product in the right-hand side of (9) converges 
uniformly on compact subsets of this line. Our first goal is to prove uniform convergence of this product 
on any compact subset of C, excluding the poles of A^(s). 

For t > and a, 2; G C satysfying Re{z) > max(0, — Re(a)), let us define 



fa{z) := In 



Tja + z) 
r(z)z'' 



9a{t) 



1 - e'"*\ 1 
1 - e-* / t' 



(10) 



It is known that for a G M, the function z G (max(0, — a), 00) 1— ?■ |/a(-2)| is completely monotone. This 
follows from [29, Theorem 4] if a > and from [28, Theorem 1] if a < 0. We also have the following 
integral representation. 



00 

fa{z) = J 5fa(t)e~^*dt, a,zeC, Re{z) > max(0, -Re(a)). 



'IV. 



which can be established using formula 8.361.5 in [12]. 
Let A be a compact subset of C. Define 

= min{Re(s) : s G A}, = max{Re(s) : s G A}. 

The convergence of the infinite product (9) is not affected by any finite number of terms. Therefore, 
without the loss of generality we can assume that 1 — ai < . If this was not true, we would take N 
large enough, so that 1 — < v~ , and investigate the convergence of the tail nn>7v(- ■ ■) of the infinite 
product in (9) . 
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For > 1 we define 

u ( \ _ rr r(/3„)r(Q„. + s-i) f 
llrK)r(/3„ + s-i) 

The function A^fc(s) is analytic and zero-free in tlie lialf-plane Re(s) > 1 — ai, and due to our assumption 
1 — «! < v~, this half-plane includes the set A. Using (11) and Lemma 1 we find that for all / > A; > 1 
and s G A, 



\HMiis))-HMkis))\ 



Y (/s-lK)-/.-l(/3n)) 



n=k+l 
oo 



oo J 
^ I 1 



-l3nt 



) dt 



n=k+l 



I 1 



?l=fc+l 

For any fixed s & A, the right-hand side in the above inequality can be made arbitrarily small if k is 
sufficiently large, since +oo as k ^ +oo. This shows that A^fc(s) forms a Cauchy sequence and, 

therefore, the infinite product in (9) converges pointwise. 

Now we need to establish the uniform convergence on the compact set A G C Using (8), we check 
that for all s in the half-plane Re(s) > 1 — ai we have M.kis) = E[X^~^], where Xk is defined by 



n=l 



Denoting v = Re(s) and using Lemma (1) and monotonicity of the function z G (max(0, — Cl), oo) I — }■ fai^Z) 
for a G M we obtain 



Mkiv) 



= exp ^^(/^„i(a„) - fv^i{f3„))^ < exp(|/„__i(ai)|). 



(12) 



The function v i— )■ /t,-i(ai) is continuous on the interval v G [v~ , v~^], therefore it is bounded on this in- 
terval. This fact combined with the inequality (12) implies uniform boundedness of all functions A^fc(s), 
/c > 1 in the vertical strip v~ < Re(s) < v~^, therefore also on the set A. By the Vitaly-Porter Theorem, 
uniform boundedness and pointwise convergence of Aik{s) imply uniform convergence of these functions. 
This proves that the product (9) converges uniformly to a function which is analytic in Re(s) > 1 — ai. 
By analytic continuation we conclude that the Mellin transform of X{a, /3) exists everywhere in this 
half-plane. □ 

In the next proposition we show that the Mellin transform of X(a, /3) satisfies two important identities, 
which will be crucial for our results on exponential functionals of meromorphic processes in section 3. 



Proposition 2. Assume that an and /3n satisfy the interlacing property (4). 
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(i) For n > 1 define an = Pn o,nd /3„ = an+i- The sequences an and Pn also satisfy the interlacing 



property (4). Define M{s) = E X{a,/3) 



\s-l 



We have the following identity 



M{s) X Mis) 



p. X .-1 > Re(s)>l-ai. 
T{ai)al 



(13) 



(a) The function Ai{s) satisfies Ai{s + 1) = (f){s)M{s) for all s E C, where is a meromorphic 
function defined as 



s-l 



nTT# 



s G C. 



(14) 



n>l 



Proof. The proof of (i) follows at once from (9). Lemma 1 implies convergence of the infinite product in 
(14), and the proof of (ii) follows easily from (9) and identity r(s + 1) = sr(s). □ 



3 Exponential functional of meromorphic processes 

Let us review some important properties of meromorphic processes, which we will need for our further 
results (see [17] for more details). The density of the Levy measure of a meromorphic process is defined 
by (3). It is clear that the Levy measure has exponentially decreasing tails. The only restriction on the 
coefficients determining the Levy measure is the convergence of the series 



n>l 



which is easily seen to be equivalent to the convergence of the integral J^x'^iT{x)dx. 

From the Levy-Khintchine formula we find that the Laplace exponent iplz) := lnE[exp(2;Xi)] of the 
meromorphic process X is given by 

tP{z) = lah^ + fiz + z^y2^^ ^ + ^2y-_^n ^ -pi < Re(z) < pi, (15) 

where a > and p G M. The function ip{z) can be analytically continued to a real meromorphic function, 
having only simple poles at points p„, — pn- The most important analytical property of meromorphic 
processes (also satisfied by hyper-exponential processes) is that for any g > the equation ipi^) = Q has 
only simple real solutions Cn? ~Cn? which satisfy the interlacing property 

••• -p2< -C2 < -Pi < -Ci < < Ci < Pi < C2 < P2 < ••■ (16) 

When we need to emphasize the dependence of these solutions on the parameter q, we will denote them 
Cn{q) and Cn{q)- Moreover, we have the following infinite product representation 

1 - ^ 1 + ^ 

q-ij{z) = ql[^—^xl[^^, zee, (17) 

n>l Pn n>l Pn 
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where the two infinite products converge due to Lemma 1. Formula (17) is the key result in the theory 
of meromorphic processes, since it allows to identify the Wiener-Hopf factors for the process X (see [17, 
Theorem 2]). 

The next theorem is the main theoretical result in this paper. We identify the distribution of the 
exponential functional Iq for a general meromorphic process. This result should be compared with [6, 
Theorem 4.3], where the distibution of the exponential functional is identified for a hyper-exponential 
process having a non-zero Gaussian component. 

Theorem 2. Assume that q > 0. Define po = and the four sequences 

C = {Cn}n>l, P = {Pn}n>l, C = {1 + Cn}n>l, P = {1 + Pn-l}n>l- 

Then we have the following identity in distribution 



where 



C = C{q):=q-'ll^^, (19) 



n>l Cn 



and the random variables X{p,Q and X{(,p) are independent and are defined by (5). The Mellin 
transform A^(s) = A^(s, q) := E[/q^^] is finite for < Re(s) < 1 + Ci and is given by 



u( ^ - TT r(Cn + i)r(Pn-i + g) L + i y r(p„)r(Cn + 1 - s) fCny~' .... 
" L\ r(Pn^i + i)r(C„, + s) [pn-i + 1 ; r(c„)r(p„ + 1 - .) [pj ' ^'^^ 



Proof. First, we note that the two pairs of sequences (C,p) and (p, C) satisfy the interlacing property 
(4), therefore the random variables X{(,p) and X(p,() are well-defined. The infinite product defining 
constanct C converges due to Lemma 1. 

Let /(s) denote the function in the right-hand side of (20). Define Mi(s) := E[(X(p, ^))*~^] and 
M2(s) := E[{X{C, p)y-^]. Using Theorem 1 we find that f{s) = C^-^Mi(s)M2(2 - s) is the Mellin 
transform of the random variable in the right-hand side of (18), and that f{s) is analytic in the strip 
< Re(s) < 1 + Ci- Formula (20) implies that /(s) is zero-free in the wider strip — < Re(s) < 1 + pi. 

We will prove that A^(s) = /(s), which would imply the identity in distribution (18). We will use the 
verification result, see [19, Proposition 2], which states that Ai{s) = f{s), provided that the following 
three conditions are satisfied: 

(i) for some 6 > 0, the function f{s) is analytic and zero free in the vertical strip < Re(s) < 1 + 6; 

(ii) the function /(s) satisfies 

/(s + l) = '-Y-Ji^)^ 0<s<9, 

where il^{s) is the Laplace exponent of the process X; 

(iii) \f{s)\^^ = o(exp(27r|Im(s)|)) as Im(s) — )■ oo, uniformly in the strip < Re(s) < 1 + 9. 
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According to the above discussion of the analytic properties of f{s), condition (i) is satisfied with 
6 = (^i. Let us verify condition (ii). Proposition 2(ii) gives us two identities 

Mi(s + l)=Mi(s)x J]^^t_l±^, M2(l-s) = M2(2-s)x s e C. 

n>l 1+Cii n>l Cn 

Combining (20) with the above identities we obtain 



fis + 1) = C-'M,is + 1)M2(1 - s) = f{s) X g-^ n I 



. + i 1 + ^ 1 - f 

n>l (n 1 + Cn 

where we have also used (19). Rearranging the terms in the above infinite product and using (17) we 
conclude that 

/(..i)^/wx..-.ni±|xi^^-i^/w, 

thus condition (ii) is also satisfied. 

Finally, let us show that condition (iii) holds true. We use Proposition 2(i) and find that 

f(s)-' = C'''—^^^^^^—xM,{s)xM2{2-s), 0<Re(s)<l + Ci, (21) 

r(s)r(Ci + 1 - s) 

where Mi(s) := E,[{X{a, l3)y~^] and M2(s) := E,[{X{ce, l3)y^^] and the sequences a, /3, a, /3 are defined 
as follows 

a„ := 1 + Cn, (3n-= 1 + Pn, ■= Pn, f^n ■= (n+1, U > 1. 

According to Theorem 1, the Mellin transform Mi{s) (resp. M2{s)) is finite in the half-plane Re(s) > — Ci 
(resp. Re(s) > 1 — pi). Since |Mj(s)| < M((Re(s)), we see that for any e > 0, the function Mi(s) (resp. 
M2(2 — s)) is uniformly bounded in the half-plane Re(s) > e — Ci (resp. Re(s) < 1 -|- pi — e). Taking 
e = ^min(Ci,pi — (i) we conclude that the function Mi(s) x M2(2 — s) is uniformly bounded in the 
vertical strip < Re(s) < 1 + Ci- 

To estimate the gamma functions in (21), we use formula 8.328.1 in [12] 

.f|y|| 



lim |r(x + i?/)|e2'^'|?/|2 = v27r, x,y eR. 

It is known that the limit exists uniformly in x on compact subsets of M (as can be easily seen from 
Stirling's asymptotic formula for the gamma function). The above formula shows that for any e > 

o(exp((7r + e)|Im(s)|)) 



|r(s)r(Ci + i-s) 



as Im(s) — )■ oo, uniformly in the strip < Re(s) < 1 + Ci. This fact combined with (21) and uniform 
boundedness of Mi(s) x M2(2 — s) shows that condition (iii) is also satisfied. Therefore, according to 
[19, Proposition 2] we have M.{s) = f{s), and this ends the proof of Theorem 2. □ 
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The results of Theorem 2 may be useful for studying self-similar Markov processes. It is well-known 
that Lamperti transformation [21] maps a positive self-similar Markov process y to a Levy process 
X, and the exponential functional of X plays a very important role in describing various properties 
of the original process Y . The last several years have witnessed a large volume of research on self- 
similar Markov process, which are constructed from stable Levy processes by conditioning or various 
path transformations (see [4, 5, 19, 20]). We note that in all of these examples (at least in dimension 
one), the Lamperti transformed process belongs to the family of meromorphic processes. Therefore, we 
hope that our results on exponential functionals of a general meromorphic process may be useful in the 
future for studying other interesting self-similar Markov processes. 

Remark 1. The results of Theorem 2 can be easily extended to the boundary case g = 0. The exponential 
functional Jq = exp(Xt)dt is well-defined if ]E[Xi] < (see [3, Theorem 1]). One can check that 
condition ?/''(0) = E[Xi] < implies Ci(q') Ci(0) = as g 0+, moreover q/Ci{q) |IE[-^i]|- Thus 
the constant C = C{q) defined by (19) converges as g — )■ 0+ to 



1 + ^ 

Pn 



n>l Cn+l(0) 



Note that Cn+i(0) > p„ > for n > 1, thus the above product is well-defined, it converges due to Lemma 
1. The random variables X(p, Q) and X(C, p) are also well-defined in the limit g — 0^, provided that we 

identify -B(i,o) = 1- 

Remark 2. The factorization of the exponential functional as a product of two independent random 
variables in (18) is in fact a very general phenomenon. It turns out that in many cases the exponential 
functional Iq has the same distribution as a product of an exponential functional of a negative of a 
subordinator and an independent exponential functional of a spectrally positive process, both of these 
processes being related to the Wiener- Hopf factors of the original Levy process X, see [24, 26]. Also, we 
would like to point out that recently Patie and Savov [25] have obtained some very strong and general 
results on the Mellin transform of the exponential functional. In particular, they show that the Mellin 
transform can be obtained as a generalized Weierstrass product in terms of the Wiener- Hopf factors of the 
process X. Since the Wiener-Hopf factors of the meromorphic process are known to be infinite products 
of linear factors (see formula (17) and [17, Theorem 2]), these results could lead to an alternative proof 
of Theorem 2. The Mellin transform of Iq can be expressed as a double infinite product, interchanging 
the order in this product and using Weierstrass infinite product representation for the gamma function 
one could obtain formula (20), probably at the cost of considerable technical details. 

4 Numerical examples 

For our numerical examples we will consider a theta process, defined by the Laplace exponent 

Lvr(V(ai-2)//3i)'' 'coth (ti ^ {ai - z) / (5i^ (22) 

+ C27r (^A/(a7+^)7/32) coth (^tta/ (^2 + z)/(5'^ 



2 

i){z) = + /iz + 7 + (-1)^' 



Ci7 



(23) 



These processes were introduced in [15] as examples of meromorphic processes for which the Laplace 
exponent can be computed in a particularly simple form. The parameter j can take two values j G {1,2}, 
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the parameter 7 is always chosen so that ip{0) = 0. Other restrictions on the parameters are 

a > 0, /i G M, Q > 0, ai > 0, A > 0. 

As was shown in [15], the density of the Levy measure of a theta process is given by (3), with parameters 
a„, dn, Pn, Pn defined as follows 

Pn = ai + Pin"^, pn = a2 + ^2^^, anpn = 2CiPin^\ dnpn = 202^2^^^ pn, u > 1. (24) 

For our numerical experiments, we will consider the following two parameter sets 

Parameter set 1: j = 1, cr = 0.1, (25) 
Parameter set 2: j = 2, a = 0.0, 

with the remaining parameters being fixed at 

p = 0.1, ci = 0.15, C2 = 0.3, «! = 02 = 1.5, (3i = (32 = 2. 

The parameter set 1 defines a process with a non-zero Gaussian component and jumps of infinite activity 
but finite variation, while the parameter set 2 defines a process with zero Gaussian component and jumps 
of infinite variation (see [15, Proposition 4]). 



4.1 Approximating the Mellin transform of the exponential functional 

First we will discuss the problem of computing the Mellin transform and the density of the exponential 
functional. This algorithm will be useful later for our discussion on pricing of Asian options. We assume 
that g > and denote the density of Ig by p{x). We start by expressing p{x) as an inverse Mellin 
transform of M{s) = IE[/^"^] 

p(x) = — [ 7W(c + iM)e~^"'"(")dM, (26) 

2-71 J 

where c can be any number in the interval (0, 1 + Ci)- We see that there are two main issues in computing 
p{x). In order to compute the Fourier integral in the right-hand side of (26), we will use Filon's method. 
This algorithm is well-known and we will refer the reader to [10, 11, 13]. However, before we can apply 
Filon's method, we need to be able to compute Ai{s) (given by an infinite product (20)) to a reasonably 
high precision. 

A naive way of computing the Mellin transform Ai{s) would be to simply truncate the infinite 
products in (19) and (20) at n = iV. After rearranging the terms, this would give us an approximation 

f ^ us-i TT r(pn-i + s) r(c„ + 1 - s) 

a. x„„ xn rK„ + ,) r(p„ + l-.) <2^) 



where 



, ^ + Pn YJ CnCn 

q ^,Pnpn 

n=l 
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and Oat is a normalizing constant, uniquely determined by the condition A^Ar(l) = 1. According to 
Theorem 2, A^iv(s) Ai{s) a.s N +00. The problem is that the convergence could be quite slow, 
and we need to find a way to accelerate it. 

Let us define the remainder term Rn{s) = A4(s)/A^Ar(s). It is clear hat Rn{s) is the tail of the 
infinite product in (20). Instead of replacing Rn{s) by 1, let us try to do a better job of approximating 
this function. 

Using Theorems 1 and 2 one can see that Rn{s) is the Mellin transform of some random variable, 
which we will denote by e*^^-*. This random variable can be obtained by taking the tail in the infinite 
products (5) and (19), however the exact form of e'-^-' is not important. What is important is that we can 
compute exactly the moments := E[(e*-^'')'^]. Note that the function Rn{s) = E[(e'-^'')*'~^] is analytic 
in the strip — < Re(s) < 1 + Ca^+i^ therefore the moments are finite for all /c < 1 + Ca^+i- Using 
the functional equation A^(s + 1) = sA^(s)/(g — ip{s)), we find 

and these quantities can be readily computed (in case if one of ip{j) is equal to zero or infinity, one can 
still compute via L'Hospital's rule). 

Our plan is to replace e^^-* (whose distribution we can not compute explicitly) by a simple random 
variable ^, with known distribution and Mellin transform, so that the first two moments of ^ match the 
first two moments of e'-^-'. We will take ^ to be a beta random variable of the second kind, which is 
defined by its density 

P(e e dx) = ^P^y--\1 + y)-'^-'dy, y > 0. 
L [a + 0) 

where the parameters a and b must be positive. The Mellin transform of ^ is given by 

^r,.-ii ^ r{a + s-i)r{b + i-s) 

One can check that if we define the parameters 

mi +1712 , 1711+1712 

a = mi 0=1H ^, (29) 

m2 — mf m2 — mf 

then we have E,[^] = mi and ]E[^^] = m2. 

Let us summarize the algorithm for approximating the Mellin transform Ai{s). We set to be a 
large number (large enough so that the condition Ctv+i > 1 is satisfied). Then mi and m2 are finite, 
and we compute these numbers using formula (28). We evaluate the parameters a and b via (29) and 
approximate the Mellin transform Ai{s) by 

. X . . / ,T(a + s - l)r(6 + 1 -s) 

We emphasize, that the right-hand side of (30) is the Mellin transform of a random variable Iq^\ which 
converges to Iq in distribution as N ^ +00, and which has the same first two moments as Ig if the latter 
exist, and "analytically continued" moments if the classical moments do not exist. 
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Figure 1: (a) The density of the exponential functional Ji with = 400 (the benchmark), (b) The error 
with iV = 20 (no correction), (c) The error with = 20 (with correction term). Solid line (resp. circles) 
represent parameter set I (resp. II). 



We illustrate the efficiency of this approximation by a numerical example. We compute the density 
of the exponential functional Ji for the two parameter sets (25), using approximation (27) with = 400. 
The graphs are shown on figure 1(a), and we will take these results as our benchmark. Then we calculate 
the density (using the same approximation (27)) with = 20, the error between these results and our 
benchmark is shown on figure 1(b). We see that the maximum error is of the order l.Oe-3. Finally, we 
perform the same calculation with A^ = 20, but now we use the approximation for the Mellin transform 
with the correction term (30). On figure 1(c) we see that the maximum error is of the order l.Oe-6, 
therefore our approximation (30) decreases the error by a factor of 1000, and it seems to be very efficient. 



4.2 Computing the price of an Asian option 

We recall that the stock price process is defined as a geometric Levy process St = SqC^* (where X is 
a Levy process started from zero and 5*0 is the initial stock price). We are working under risk-neutral 
measure, so that the process Stexp{—rt) is a martingale under measure P). This condition can be 
expressed in terms of the Laplace exponent of the process X by the equation ipi^) = (note that a 
condition pi > 1 is necessary to ensure ip{l) < +00). In order to satisfy this risk-neutral condition, we 
will choose accordingly the value of /i (which is responsible for the linear drift of the process) in formula 
(22). The price of a fixed strike Asian option is given by the expectation (1). Let us introduce the 
function 



f{k,t) :=E 



Since C{So, K,T) = exp(— rT) x Sq x f{K/SQ,T), we can study an equivalent problem of computing 
numerically the values of the function f{k,t). 

Our first goal is to connect the function f{k, t) with the Mellin transform of the exponential functional 
Iq. This result was given in [6, Theorem 2.1], however we reproduce the main steps for the sake of 
completeness. We see that for any q> r 



h{k,q) 



POO 

:=q e"''fik,t)dt = E[iIq-k)+] 
Jo 



(31) 
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where e{q) is an exponential random variable with mean 1/q (independent of X). Note that h{k, q) < +00 
for q > r, since E [(Jg — k)~^~^ < E[/g] = {q — r)~^. Next, let us check that for g > r we have (i > 1. This 
is clear, since Ci is the smallest positive solution to the equation tp{z) = q, and we have the risk-neutral 
condition = r. Since ip{z) is convex and ip{0) = 0, we find that > 0, thus the solution to 

equation ip{z) = q must be greater than 1. Having established that (i > 1, we verify that for q > r and 
< s < ^1 — 1 we have 



h{k,q)k'~Uk 



E 
E 



" POO 

/ {Ig - kf k'-^dk 

Jo 



{Ig - k) k'-'dk 



(32) 



E[jfi] M{s + 2,q) 



s(s + l] 



s(s + l] 



where in the first step we have used Fubini's theorem. Note that the right-hand side is finite for < s < 
(i — 1, since ^A{s) is finite in the strip Re(s) G (0, 1 + Ci) (see Theorem 2). The two formulas (31) and 
(32) are the foundations for our algorithm to compute the price of an Asian option. 

Algorithm 1: approximating the Mellin transform of the exponential functional. 

This algorithm is based on inverting Laplace and Mellin transform in (31) and (32) and approximating 
Ai{s) by the algorithm presented in section 4.1. First, for d2> r and g = ^2 + we compute h{k, q) as 
the inverse Mellin transfom 



h{k,q) 



k-'^^ r M{di + iv + 2,q) 
27T J {di + iv){di + iv + 1) 



(33) 



where di G (0,Ci(c?2) — !)• Second, we compute f{k,t) as the inverse Laplace transform, which can be 
rewritten as the cosine transform 



fik,t) 



Re 



h{k, d2 + iu) 
d2 + iu 



cos{ut)du. 



(34) 



We set di = d2 = 0.25 and truncate the integral in (33) (resp. (34)) so that the domain of integration 
is —100 < V < 100 (resp. < m < 200), and use Filon's method [10, 11] with 400 discretization points to 
evaluate each of these integrals. The Mellin transform is computed using the approximation algorithm 
presented in section 4.1, with the Mellin transform truncated at N terms (we will set G {10, 20, 40, 80} 
in our computations). Computing the Mellin transform requires computing 2A^ solutions to the equation 
ip{z) = q. An efficient algorithm of how to find these solutions for both real and complex values of q was 
presented in [14]. 

We would like to point out that our algorithm requires computing the values of Ai{s, q) for complex 
values of g, while Theorem 2 was established only for g > 0. To extend the results of Theorem 2 for 
complex values of g, the main step would be to establish the uniform convergence of the infinite product 
in the right-hand side of (20) for complex g, for which we will need some additional information on the 
behavior of the solutions Cn, — Cn to the equation 'ip{z) = q for complex g. A rigorous discussion of this 
question is beyond the scope of the present article, and we will leave it to future work. 

Algorithm 2: approximating theta-process by a hyper-exponential process. 

Our second algorithm is based on approximating the theta-process X by a hyper-exponential process 
X = X^^\ for which the Mellin transform of the exponential functional can be computed explicitly. 
There is a natural way of obtaining such an approximation for any meromorphic process: We can simply 
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truncate both infinite series defining the density of the Levy measure in (3) at terms (see [8] for 
another approximation technique). This procedure gives us a Levy process X with hyper-exponential 
jumps, whose Laplace exponent is given by 



N 

- ~2 2 , ~ I 2 

-a z + HZ + z 



1 Pn{Pn - Z) 
n=l ^ ' 



+ z- 



N 



(35) 



see (15) for comparison. In our case, the coefficients a„, p„, a„ and p„ are defined by equation (24). 
We will choose a so that the variance of Xt matches the variance of Xj, which is equivalent to requiring 
'?/'"(0) = V'"(0). The parameter fx is then specified by enforcing the risk-neutral condition = r. 

Following this simple procedure we obtain a sequence of hyper-exponential processes X^^\ which will 
converge in Skorohod's topology to the process X. Now we can compute the price of the Asian option, 
with the driving process X, following the same procedure as in algorithm 1. The only difference is that 
the Mellin transform of the exponential functional Iq = J^^"^^ exp(Xj)dt is computed as 



M{s,q) = E 



rs-l 



a X 



~o\ 1— s 



N N+l 

nr(p, + s) n r(i + o-5) 



X r(s) X 



X 



N+l N 

n r(0 + s) n r(i + Pj - s) 



where a = a{q) is chosen so that Ai{l,q) = 1. This expression for the Mellin transform was obtained in 
[6]. Here (n and — Cn are solutions to the equation ipi^) = and since ipi^) is a rational function, it is 
easy to see that there will be exactly + 1 positive and A^ + 1 negative solutions to ipi^z) = q (see [6] 
for all details). 

Algorithm 3: Monte-Carlo simulation. 

We will also check the accuracy of the previous two algorithms by computing the price of the Asian 
option by a simple Monte-Carlo simulation. We approximate the theta-process X = {Xt}o<t<T by a 



random walk Z = {Zn}o<n<Aoo with Zq = and the increment Z, 
Asian option is approximated then by the following expectation 



n+l 



X' 



r/400- 



The price of the 



^ 400 



400 



K 



n=l 



which we estimate by sampling 10^ paths of the random walk. In order to sample from the distribution 
oiY := Zn+i — Zni we compute its density Py{x) via the inverse Fourier transform 



27r 



where E [e'""^] = E [e'^^^^/^ooj = exp ((T/400)'?/'(i2;)) and the Laplace exponent iIj{z) is given by (22). 
Again, in order to compute the inverse Fourier transform, we use Filon's method with 10^ discretization 
points. 

We compute the price of the Asian option with the initial stock price 5*0 = 100, interest rate r = 0.03, 
maturity T = 1 and strike price K = 105. We consider the two theta-processes defined by parameter sets 
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Algorithm 1, price 


Algorithm 1, time (sec.) 


Algorithm 2, price 


Algorithm 2, time (sec.) 


10 


4.724627 


1.6 


4.720675 


1.2 


20 


4.727780 


2.8 


4.728032 


1.8 


40 


4.728013 


4.8 


4.728031 


3.4 


80 


4.728029 


9.2 


4.728031 


7.1 



Table 1: The price of the Asian option, parameter set I. The Monte-Carlo estimate of the price is 4.7386 
with the standard deviation 0.0172. The exact price is 4.72802ihl.0e-5. 





Algorithm 1, price 


Algorithm 1, time (sec.) 


Algorithm 2, price 


Algorithm 2, time (sec.) 


10 


10.620243 


1.6 


10.621039 


1.2 


20 


10.620049 


3.0 


10.620171 


2.2 


40 


10.620037 


4.8 


10.620054 


3.6 


80 


10.620036 


9.6 


10.620039 


7.4 



Table 2: The price of the Asian option, parameter set II. The Monte-Carlo estimate of the price is 
10.6136 with the standard deviation 0.0251. The exact price is 10.62003±1.0e-5. 



I and II (see (25)). Note that the parameter /x is not equal to 0.1 anymore, as it has to be determined 
by the risk-neutral condition = r. The results of our computations are presented in tables 1 and 2. 
The code was written in Fortran90, and all computations were performed on a basic laptop (CPU: Intel 
Core i5-2540M 2.60GHz). 

The results presented in tables 1 and 2 show that both algorithm 1 and algorithm 2 perform very 
well, and seem to converge quickly to the true value of the option. The exact price was computed with 
N = 160 and 1600 discretization points for the two integrals in (33), (34). By experimenting with other 
values of the above parameters, we are convinced that these values are correct to withing l.Oe-5. Both 
of these algorithms are very efficient; the CPU time seems to be comparable with the results of Cai 
and Kou on hyper-exponential processes (see [6]). Note that there is a substantial difference between 
algorithm 1 and algorithm 2, as one is based on approximating the Mellin transform of the exponential 
functional, and the other is based on approximating the underlying Levy process by a hyper-exponential 
process, for which the Mellin transform of the exponential functional can be computed explicitly. Yet the 
results produced by both of these algorithms agree up to five decimal digits, which is a good indicator 
that these results are indeed correct. The algorithm based on Monte-Carlo computation also produces 
consistent results, however these estimates are much less accurate and require CPU time on the order of 
several minutes. 

We have performed similar numerical experiments for many other values of parameters of the un- 
derlying theta-process, as well as for different maturities and different strike prices. Qualitatively, the 
results seem to be consistent with the ones presented in tables 1 and 2. Algorithm 2 is very efficient for 
the parameter set I, and gives high accuracy even for relatively small values of N. This should not be 
surprising, as in this case we have a theta-process with jumps of finite variation (but infinite activity), 
and it is intuitively clear that processes with compound Poisson jumps can provide a good approximation 
for such processes. On the other hand, parameter set II corresponds to a theta-process with jumps of 
infinite variation and zero Gaussian component, and here our approximation by a compound Poisson 
process with a non-zero Gaussian component is bound to be less precise. Nevertheless, algorithm 2 works 
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quite well in all cases, and, in our opinion, it may be a preferred method to compute prices of Asian 
options for any meromorphic processes (though testing this algorithm on other meromorphic processes, 
such as beta-processes [14], would be worthwhile). 

Algorithm 1 has comparable performance (though it is a little slower than algorithm 2). It has one 
potential advantage compared with algorithm 2. Imagine that we want to compute the price of the Asian 
option for = 20 and then to check whether we have sufficient accuracy by doing the same computation 
with A^ = 40. Algorithm 2 would require re-computing everything, since the Laplace exponent ip of the 
hyper-exponential process will change, thus all the numbers (n, —(n (the solutions to equation ip{z) = q) 
will be different. This is not the case for algorithm 1: here the Laplace exponent ip{z) does not depend 
on A^, thus the roots Cn and with 1 < n < 20 will not be affected. Therefore, we only need to compute 
the new roots for 21 < < 40, which will need fewer computations. The same idea can be applied for 
the evaluation of the Mellin transform, where the results for A^ = 20 in (27) can be stored in memory, 
and only the remaining finite product of gamma functions with 21 < n < 40 have to be evaluated. 
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